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ABSTRACT 

The propagation of charged particles, including cosmic rays, in a partially ordered magnetic 
field is characterized by a diffusion tensor whose components depend on the particle’s Larmor 
radius i?L and the degree of order in the magnetic field. Most studies of the particle diffusion 
presuppose a scale separation between the mean and random magnetic fields (e.g., there be¬ 
ing a pronounced minimum in the magnetic power spectrum at intermediate scales). Scale 
separation is often a good approximation in laboratory plasmas, but not in most astrophysical 
environments such as the interstellar medium (ISM). Modern simulations of the ISM have 
numerical resolution of order 1 pc, so the Larmor radius of the cosmic rays that dominate in 
energy density is at least 10® times smaller than the resolved scales. Large-scale simulations 
of cosmic ray propagation in the ISM thus rely on oversimplified forms of the diffusion tensor. 
We take the first steps towards a more realistic description of cosmic ray diffusion for such 
simulations, obtaining direct estimates of the diffusion tensor from test particle simulations 
in random magnetic fields (with the Larmor radius scale being fully resolved), for a range of 
particle energies corresponding to 10’ ^ Rh/lc ^ 10®, where Ic is the magnetic correlation 
length. We obtain explicit expressions for the cosmic ray diffusion tensor for R^/lc ^ 1, 
that might be used in a sub-grid model of cosmic ray diffusion. The diffusion coefficients 
obtained are closely connected with existing transport theories that include the random walk 
of magnetic lines. 

Key words: magnetic fields - diffusion - MHD - turbulence - methods; numerical - cosmic 
rays 


1 INTRODUCTION 

The theory of propagation and confinement of cosmic rays is 
largely based on analytical calculations (mostly quasilinear, ap¬ 
plicable if magnetic fluctuations are much weaker than the back¬ 
ground magnetic field) of the scattering of charged relativistic 
particles by magnetic inhomogeneities (Berezinskii et al. 1990; 
Schlickeiser 2002; Kulsrud 2004; Strong et al. 2007). Analytical 
results have been verified and extended by simulations of the mo¬ 
tion of test particles (Giacalone & Jokipii 1994, 1999; Casse et al. 
2001; Candia & Roulet 2004; Plotnikov et al. 201 1) in model mag¬ 
netic fields B represented as the sum of a mean field Bo and fluc¬ 
tuations b, 

B^Bo + b. (1) 

The main goal of the test particle simulations has been the verifi¬ 
cation of the analytical results and their extension, especially to the 
case of finite magnetic fluctuations (e.g., Shalchi & Dosch 2009), 
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often focusing on the dependence of the cosmic ray propagation 
parameters on the particle energy. 

The results are formulated in terms of the cosmic ray diffusion 
tensor 

Kij = K^Sij + (A'li - K^)BiBy , (2) 

where B = B/\B\ is the unit vector in the magnetic field direc¬ 
tion, and it'll and K± are, respectively, the cosmic ray diffusivities 
along and across the mean magnetic field (e.g., Giacalone & Jokipii 
1999, and references therein). Extensive numerical results have 
been obtained for a wide range of the particle rigidity, 10“^ < 
P < 10®, where 

p = Rl/Ic ~ E/q, 

with 

iin«3.3xl0®^cm(^^) sinP (3) 

the Larmor radius of the particle’s gyration in the magnetic field 
(here given for ultra-relativistic protons), Ic the correlation length 
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of magnetic fluctuations, E and q the total particle energy and elec¬ 
tric charge, and 0 the angle between the magnetic field and the par¬ 
ticle velocity (the pitch angle). The results of test particle simula¬ 
tions are in approximate agreement with theoretical estimates and 
can be fitted with the form (Section 2.3 of Strong et al. 2007) 


if|l 


KO 


E ygg 

IGeVj ’ 


(4) 


with an accompanying expression for iT||/iT^, where bo is the 
root-mean-square strength of magnetic fluctuations, n = 1/3-1/2 
depending on the magnetic power spectrum, and kq — (3-5) x 
10^®cm^s“^ for the Galactic cosmic rays near the Sun. The 
spectrum of magnetic fluctuations in the ISM extends from the 
outer scale of order 100 pc ^apparently down to at least 10® cm 
(Armstrong et al. 1995). Therefore, the Larmor radius of cosmic 
ray particles of GeV energies is well within the range of magnetic 
turbulent scales. 

The diffusion of charged particles in a random magnetic field 
is a result of two dominant effects. Particles gyrating in a magnetic 
field are scattered by magnetic fluctuations — which may be part 
of the externally driven magnetic spectrum or be self-generated by 
the streaming instability (Kulsrud 2004) — at a scale comparable 
to i?L. Due to the scattering, the particle propagation both along 
and across the magnetic field become diffusive. The divergence of 
the field lines of a random magnetic field that guides the particles 
further enhances their diffusion. As a result, the propagation of cos¬ 
mic rays is sensitive to both local and global properties of the mag¬ 
netic field. Sonsrettee et al. (2015a) present a detailed analysis of 
the random walk of magnetic lines with vanishing mean magnetic 
field (see also Snodin et al. 2013; Sonsrettee et al. 2015b). 

The dominant contribution to the scattering of charged parti¬ 
cles is provided by resonant magnetic fluctuations at a scale close 
to the particle’s Larmor radius given by Equation (3). Magnetic 
fluctuations at such a small scale can be measured or estimated 
only in rare cases, e.g., in situ in the interplanetary space. This 
makes it difficult to apply the theory of cosmic ray diffusion to 
interstellar or intergalactic cosmic rays at their dominant energies, 
E ~ (1-10) X 10® GeV per nucleon, since little beyond crude 
statistical properties of random magnetic fields at their outer (or 
energy-range) scale is known with any degree of confidence. This 
scale is of the order of 3 x 10^® cm (100 pc) in the interstellar space 
and 3 x 10^^ cm (10 kpc) in galaxy clusters, by far larger than the 
resonant scale of the cosmic ray particles. 

More importantly, models of cosmic ray propagation rely on 
the representation (1) for magnetic field, which is meaningful and 
useful only in the case of a clear scale separation between the mean 
magnetic field and its fluctuations. Under such a scale separation, 
the range of scales (spatial and/or temporal) over which the mean 
magnetic field varies is clearly distinct from that of the random 


^ This estimate of the largest scale in the interstellar inertial range follows 
from analyses of the interstellar velocity field (page 148 in Ruzmaikin et al. 
1988, and references therein). Observations of interstellar turbulence in var¬ 
ious tracers, e.g., Faraday rotation (Minter & Spangler 1996) or fluctua¬ 
tions of synchrotron intensity (Haverkorn et al. 2008; lacobelli et al. 2013), 
which depend on products of fluctuating quantities, unsurprisingly sug¬ 
gest smaller values of the turbulent scale as long as they are not inter¬ 
preted in terms of the statistical properties of each of the physical variables 
involved (thermal and cosmic-ray electron densities, magnetic field, etc.) 
(Stepanov et al. 2014). Such an inteipretation is hampered by the fact that 
the variables can be correlated or anticoiTelated to a poorly known degree 
(Beck et al. 2003; Stepanov et al. 2014). 


magnetic field; we note that the mean magnetic field does not need 
to be uniform (Gent et al. 2013; Farge & Schneider 2015). For ex¬ 
ample, scale separation can manifest itself as a local minimum in 
the magnetic power spectrum at a scale larger than the outer turbu¬ 
lent scale (see Gent et al. 2013, for a discussion). Clear separation 
of the mean magnetic field and fluctuations is typical of laboratory 
plasmas and the solar wind (Goldstein et al. 1995; Howes 2015), 
where the mean and fluctuating magnetic fields vary in space and 
time at widely separated scales. However, there are neither theo¬ 
retical nor observational reasons to expect that the spectra of inter¬ 
stellar or intergalactic magnetic fields would have such a feature. 
In any event, the existence of such a separation at a scale compa¬ 
rable to a certain value of i?L would justify Equation (1) only over 
a relatively narrow range of particle energies, rather than for the 
wide range of i?L required in a propagation model. A consistent in¬ 
terpretation of Equation (4) could be to use it as a scale-dependent 
diffusivity with the scale-dependent ‘mean’ magnetic field defined 
as Bq (k) = (k') Ak'/k' or similarly, but then the propagation 

of cosmic rays has to be described with an integral equation which 
can be reduced, under certain conditions, to a partial differential 
equation of a more complicated form than the diffusion equation 
(e.g.. Section 5.4 in Bakunin 2008, see also Snodin et al. (2006)). 

As a result, large-scale simulations of cosmic ray propagation 
in the galaxies use the forms (2) and (4), or even a scalar diffu¬ 
sivity with kq simply taken to be constant or a heuristic function 
of position designed to model galactic spiral arms, disc and halo. 
This is true even of such advanced cosmic ray modelling tools 
as the GALPROP (Strong & Moskalenko 2001; Strong et al. 2007; 
Vladimirov et al. 2011), COSMOCR (Miniati 2001) and DRAGON 
(Evoli et al. 2008; Gaggero et al. 2013; Maccione 2013) codes; see 
also Miniati (2007); Hanasz et al. (2009); Siejkowski et al. (2014). 
This is an oversimplification since both parameters that control cos¬ 
mic ray diffusion, the ratio of the Larmor radius to the turbulent 
correlation scale and the strength of magnetic fluctuations, depend 
on all MHD variables including gas density, velocity and tempera¬ 
ture. As a result, the relation of the local energy density of cosmic 
rays to the ISM parameters, a question of primary importance in 
radio astronomy, remains uncertain (e.g., Stepanov et al. 2014, and 
references therein). 

One of the goals of this paper is to develop approaches to a 
more physically detailed modelling of cosmic ray propagation in 
numerical magnetohydrodynamic (MHD) simulations of the inter¬ 
stellar medium (ISM). Since the Larmor radius of cosmic ray par¬ 
ticles is by far smaller than the numerical resolution of any ISM 
simulations, which tend to resolve scales down to the order of a 
parsec, a subgrid model of cosmic ray diffusion is required. Most 
test-particle simulations of cosmic ray diffusion in the context of 
the ISM employ a synthetic random magnetic field specified in con¬ 
tinuous space, as in Equation (9) below (see, however, Casse et al. 
2001; De Marco et al. 2007; Beresnyak et al. 2011), whereas the 
ISM simulations use a numerical grid with a certain spatial resolu¬ 
tion. This should not prevent such results from being used in grid- 
based simulations, but we consider test particles motion in both 
continuous and discrete spatial domains in Section 2.1 to appreci¬ 
ate their limitations and select optimal parameter ranges. 

Analytical studies of the propagation of charged particles 
in random magnetic fields have produced deep insights (see 
Casse et al. 2001, for a useful review), many of which show con¬ 
sistency with test particle simulations. However, such analytical re¬ 
sults are not particularly accurate over a wide range of parameters, 
presumably due to the assumptions involved in their derivation. 
Here we adopt an empirical approach and aim at deriving approx- 
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imate asymptotic forms of the diffusion tensor from test particle 
simulations, which may be useful in a wide range of applications. 
The studies of cosmic ray propagation with test particles often fo¬ 
cus on energetic particles whose Larmor radius is comparable to or 
exceeds the correlation length of magnetic fluctuations. However, a 
sub-grid model of cosmic ray diffusion requires the opposite limit, 
of a Larmor radius much smaller than the scales at which magnetic 
field is known. This limiting case is our main concern in this pa¬ 
per: we attempt to develop expressions for the cosmic ray diffusion 
tensor in a partially ordered, random magnetic field that can be ex¬ 
trapolated to scales much larger than the particles’ Larmor radius. 
Thus, in order to deduce the limiting behaviour, we perform test 
particle simulations in magnetic fields containing a large range of 
scales, but where the particle Larmor radius is well resolved nu¬ 
merically and is smaller than the correlation length of the magnetic 
field. The obtained components of the cosmic ray diffusion tensor 
can then be applied to much larger scales, depending on the particle 
momentum and the magnetic field strength on those scales. 


2 TEST PARTICLE SIMULATIONS 


Most test particle simulations employ a synthetic random magnetic 
field, either composed of a superposition of static plane waves or 
specified on a discrete mesh. Earlier simulations have demonstrated 
that the propagation of charged particles in a random magnetic 
field is very sensitive to subtle details of the magnetic field struc¬ 
ture, such as the density at which the wave vector space is pop¬ 
ulated with magnetic modes and the range of wave numbers in¬ 
volved, etc. (Casse et al. 2001; Parizot 2004; Fatuzzo et al. 2010; 
Plotnikov et al. 2011). Here we describe test particle simulations in 
which two quite different models of an isotropic turbulent magnetic 
field are used and carefully compared in order to better understand 
these subtleties. The magnetic fields are taken to be stationary. We 
neglect electric fields since they are expected to have a negligible 
effect on relativistic particle energy as the acceleration time is likely 
to be much greater than the diffusion time (Fatuzzo et al. 2010). 

In one model used here, the magnetic field is defined at the 
particle’s location using an explicit algebraic expression. The other 
model constructs the magnetic field on a regular mesh, as in MHD 
simulations. Equations of motion are then solved for a large num¬ 
ber of particles of a given Larmor radius, over a number of mag¬ 
netic field realisations. The elements of the diffusion tensor are then 
obtained as functions of time from the mean-squared particle dis¬ 
placements, with averaging over the particle trajectories and mag¬ 
netic field realisations. At a suitably large time, where the diffusion 
coefficients have settled to steady values, asymptotic diffusion co¬ 
efficients can be obtained. 

Our model magnetic field is of the form (1) where Bo = Bqz 
is a uniform mean field and 6 is a random component, with (b) = 0 
(where angular brackets denote suitable averaging). It is convenient 
to parameterize the turbulence level with 


V = 


fei + Bi ’ 


(5) 


where 0 < »7 < 1 and fog = we note that {bo/Bo)^ = 77/(1 — 
77 ). Here rj — 0 corresponds to a uniform, unperturbed magnetic 
field, whilst 77 = 1 implies no mean field, with the magnetic field 
being dominated by turbulent small-scale components. The random 
magnetic field is assumed to have an isotropic power spectrum 


M{k) = Mo 


0 

(fc/fco)-“ 


for k < ko, 
for k > ko, 


( 6 ) 


where Mo is a normalization constant related to the desired tur¬ 
bulent magnetic energy density &o, and ko = 27r/L, with L the 
largest (outer) scale of turbulence. Therefore, the spectral power 
of the magnetic field is localised at well-separated wave numbers, 
k > ko due to the random field, and k — 0 due to the mean (in this 
case, uniform) part. 

We consider two values of the spectral index, s = 5/3 or 
3/2, to explore its effects on the particle transport. The correlation 
length Ic of an isotropic field is given by (Monin & Yaglom 1975) 

, n!f’k-^M{k)dk 

" 2 /g°°M(fc)dfc ’ ^ 

with L = O.IL = 7 r/( 5 fco) for s = 5/3, and L = 7 r/( 6 fco) for 
s = 3/2. Our numerical implementations of the random magnetic 
field have discrete values of the wave number k and extend to some 
maximum k — fcmax- As a result, the effective value of the corre¬ 
lation length differs by up to a few per cent from the above values 
obtained for fcmax —>■ 00 . We also considered a magnetic spectrum 
of the form 

M(fc) = Afo (fc/fe6)2]-/2+2 ’ 

which is similar to Equation ( 6 ), but has the peak in the energy 
spectrum near some chosen kb rather than at ko- 


2.1 Magnetic field models 


The first approach is to represent the magnetic field as a sum of a fi¬ 
nite number of plane-wave modes with random polarizations, wave 
vectors and phases. It can be shown (e.g., Batchelor 1953) that the 
resulting field is spatially homogeneous and isotropic in the limit 
of an infinite number of modes. Giacalone & Jokipii (1994) imple¬ 
mented this idea which has subsequently been used in many stud¬ 
ies involving test particle simulations. Here we use a very similar 
method based on a time-independent version of the velocity field 
used by Fung et al. (1992) to model synthetic turbulence (and often 
used in modelling of turbulent flows by other authors). We take 


N 

b{x) = ^ [C. cos{kn ■ x) + Dn sin(fe„ ■ x)] , (9) 

n = l 


where are randomly oriented wave vectors, and Cn and Dn are 
random vectors confined to the plane perpendicular to to ensure 
that V • b = 0. With 


\Cn? = \Dn\^ = 


Ak„M{kn) 
AknM{kn) ’ 


( 10 ) 


where Akn is a suitable wave number spacing as defined below, the 
desired energy spectrum is obtained. The random angle between 
the amplitude vectors controls the random polarization of the plane 
waves. An efficient way to model a power-law spectrum is to select 
geometrically spaced wave numbers kn between ko and fcmax: 


kn — ko 



(n-l)/(N-l) 


( 11 ) 


where 1 < n < (so that fci = ko and fcjv = fcmax). In Equa¬ 
tion (10), we then take Afci = (fc 2 — fci), Afcjv = (fcjv — fciv-i), 
and otherwise Afc^ = (fcn+i — fcn-i)/2, for 2 < < A — 1. For 

large N, the resulting magnetic field is essentially non-periodic, 
which is potentially useful in detecting periodicity effects when 
compared with the second (periodic) model described below. An 
appropriate value for N depends on the context; we examine this in 
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Appendix A. This model of a random magnetic field will be hence¬ 
forth referred to as the continuum model. 

The second model implements a spatially periodic random 
magnetic field on a regular, three-dimensional Cartesian mesh, 
similarly to Casse et al. (2001). First, a regular grid in the three- 
dimensional wave number space is populated with complex Fourier 
components b{k) given by 

b{k) = -f ib2{k)d’^^^% (12) 

where bi and b 2 are random, real unit vectors perpendicular to both 
k (to ensure that V • b = 0) and each other, so that bi is chosen 
at random and b 2 — k x bi, where the hat denotes a unit vector. 
Flere 4>i{k) and (f> 2 {k) are distinct random phases that control the 
polarization of the magnetic field modes. We note that results pre¬ 
sented below, at least for 77 = 1 , are not particularly sensitive to 
the random phases, and taking (fii (fc) = 02 (k) = 0 changes the 
results very little: the randomness in the choice of bi is sufficient 
to ensure randomness in the relative phases of different modes. We 
take A{k) oc ^ M{k)/k to set the correct energy spectrum. We 
enforce Hermitian symmetry on b(fc) so that the inverse Fourier 
transform of Equation (12) yields a real, random magnetic field 
which we then normalise to the desired energy density bo- In what 
follows, this model is referred to as the discrete model. The appro¬ 
priate grid size is discussed in Appendix A where the two magnetic 
field models are compared. We select a sufficiently small grid spac¬ 
ing in the discrete model, and a sufficiently large maximum wave 
number fcmax in the continuum model, to ensure that the magnetic 
field is fully resolved at the Larmor radius scale. 

There are a few important subtleties in the numerical imple¬ 
mentation of random magnetic fields on a mesh (i.e. in the dis¬ 
crete model), arising from the requirements of isotropy and ran¬ 
domness of phases (to a given accuracy). To ensure these at the 
outer scale fco, there should be a sufficiently large number of modes 
at this wave number. In other words, there should be a sufficiently 
large sphere of (approximately) this radius in fc-space. Failing to 
do this may affect the particle diffusion significantly. Therefore, 
the size of the computational domain in fc-space should be larger 
than fcmax/fco, and, correspondingly, the domain in physical space 
should be larger than L = 27r/fco. If the isotropy and the random¬ 
ness of phases have been achieved at the outer scale, smaller scales 
do not represent a problem in this respect. In addition to magnetic 
modes at the resonant wave number k « 2 tt/Ri^, previous studies 
have shown that the density of nearby wave numbers can be impor¬ 
tant for effective particle scattering with some types of magnetic 
field (e.g., Mace et al. 2012). Similar consideration may be impor¬ 
tant for both of our magnetic field constructions. We carefully test 
our numerical implementations of the magnetic field for isotropy 
at all relevant scales including 27r/fco to ensure that the modes of 
a given wave number have randomly distributed phases. It seems 
that taking a mesh size Nx = Ny = Nz > 4fcmax/fco is sufficient 
if fcmax is set to be the radius of the largest sphere that fits within 
the wave number domain. Larger mesh sizes may be preferable 
when i?L ^ L, but in that case one might compensate by reduc¬ 
ing fcmax/fco, since the small scales should have little influence on 
the particle transport (Plotnikov et al. 2011). 

We use a trilinear interpolation from the nearest grid points to 
obtain the magnetic field between them. On one server with 256 GB 
of RAM we could achieve a resolution of up to 2048® mesh points. 
Larger resolutions, involving a mesh distributed in memory across 
multiple compute nodes, are much more computationally expensive 
and are not likely to provide any additional benefit in the present 
study. A magnetic field realization can be computed in much less 


than a minute of CPU time at a resolution 1024®. It is desirable to 
use a large number of the magnetic field realisations, with as few 
particles per realisation as possible. For example, at 1024® around 
30 realizations, each with about 50-100 particles can be performed 
within a reasonable time. However, in such a time, more particles 
per realisation have to be used as the mesh size increases in order to 
obtain sufficient statistics. Once the magnetic field has been com¬ 
puted, solving the particle equations of motion is relatively quick. 

When magnetic field is defined using Equation (9) (i.e. in the 
continuum model), very little computer memory is needed because 
the field is computed only at the location of a particle. In this case, 
the number of modes N is limited by the computer time available 
because the field is computed repeatedly at each time step when the 
equations of motion are solved. 


2.2 Equations of test particle motion 


The position x of an energetic charged particle in magnetic field B 
is governed by the Newton-Lorentz equations. 


dx dv q 

——= V, —— = —^v X B 

dt dt jmc 


(13) 


where v is the particle velocity, q and m are the particle charge and 
rest mass, 7 = \/l — jc? is the Lorentz factor, v = |u|, and 
c is the speed of light. Since the electric field has been neglected, 
particle energy is conserved. It is useful to define dimensionless 
variables 






t 

to 


with to = L/vo and Bo = {Bq + 60 )in terms of L and an 
appropriate reference speed uq; for relativistic particles, uq ~ c is 
an appropriate choice. The resulting dimensionless form of Equa¬ 
tions (13) is given by 



dv' 

dt' 


= av' X B' , 


(14) 


where 


qLBo 

"fmcvo 


(15) 


In this work we examine how the particle transport depends on the 
characteristic Larmor radius. 


_'ymvc_ pc _pcy/p 
\q\Bo \q\Bo k| 60 ’ 

where p = 'ymv is the particle momentum. We note for the sake 
of the reader’s convenience that i?L is related to the particle kinetic 
energy via Ek = sj{R'LqBoY + vrBcA — mc^, or Ek oc 7 ?l for 
relativistic protons in a microgauss-strength magnetic field. The di¬ 
mensionless Larmor radius is given by 


R'l — RcjL — v' I a . 


We note that the Larmor radius is defined here with respect to the 
total rather than mean magnetic field strength. 

Individual particles have random initial positions within a 
cube of approximate length L, the outer scale of the turbulence. The 
particle velocity is given a random initial direction with |t;'| = 1 . 
We then solve equations of motion numerically for Hq = 1 and set 
by selecting a as desired. 

We solve Equations (14) numerically, primarily using the fifth- 
order method of Cash & Karp (1990), which allows for an adaptive 
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step size via the error estimate of the fourth-order solution. Addi¬ 
tionally, in order to check if precise energy conservation is essen¬ 
tial, we have sometimes used the highly conservative method of 
Boris (1970), which is widely used in particle-in-cell (PIC) plasma 
simulations. We find no appreciable difference in results between 
the two methods; however, the former is more efficient. Over all 
results reported here, the maximum particle kinetic energy change 
over the diffusion time fmax - which is typically ( 10 ^- 10 ^ )fo. as is 
seen in Figure A1 - is around 0.1 per cent (but is often much less). 


2.3 Diffusion coefficients 


Time-dependent diffusion coefficients were obtained from test par¬ 
ticle simulations using 


K.xa:{t) = = (Vx:{t)Ax) , (17) 

where Ax is the particle displacement in the i-direction over 
some time t, and Vx (t) is the i-component of the particle veloc¬ 
ity at the end of the displacement. The angular brackets in Equa¬ 
tion (17) denote averaging over displacements for multiple time 
intervals of length f on a specific particle trajectory, over many 
particles and several (ideally, equally many) magnetic field re¬ 
alisations. Similar expressions were used to evaluate Kyy{t) and 
Kzz(t). In this work we are mainly interested in the asymptotic 
limit, Kij = limt_»oo which can be approximated by tak¬ 
ing at a suitably large time. In this limit, Equation (17) is 

equivalent to 


f^XX 


lim 

t—FOO 


2t 


(18) 


which has been employed in other works (Section 1.3.1 of Shalchi 
2009, and references therein). As illustrated and discussed in Ap¬ 
pendix A, the two prescriptions are equivalent but the former is 
more convenient in simulations, especially when the time depen¬ 
dence of the diffusivity is as smooth and simple as in our case. In 
particular, the asymptotic diffusion coefficients are obtained from 
Equation (17) at an earlier computational time than with Equa¬ 
tion (18). 

Here we have no real interest in the off-diagonal elements Hij , 
where i 7 ^ j, since they describe systematic drifts of cosmic-ray 
particles that are not likely to affect their large-scale diffusion (Sec¬ 
tion 1.6 of Shalchi 2009). We use the difference between the nu¬ 
merically obtained K^xf) and Kyy{t) as an estimate of the error in 
the diffusion coefficients; typically we obtain sufficient statistics to 
have the maximum running difference of less than three per cent. 
Note that for the diagonal elements of Kij — but not, in general, 
for the off-diagonal ones (Shalchi 2011) — the above relations pro¬ 
duce the same result as the standard Taylor-Green-Kubo formula 
(Taylor 1922; Green 1951; Kubo 1957) 

^00 

= / {vi{0)vj{t))dt, (19) 

Jo 

where Vi{t) is the velocity component in the direction Xi at time t. 

The diffusion coefficients Kij thus obtained are related, via 
a relation similar to Equation (2), to the diffusivities k,± and Ky 
across and along the magnetic field at the largest scale available, 
which is comparable to L in our implementations of the random 
magnetic field. In contrast, the parallel and perpendicular diffusiv¬ 
ities 7T|| and K± of Equation (2) are defined with respect to Bq, 
the magnetic field averaged over scales larger than the turbulent 
ones. Therefore, we call Kij, Ky and k± the local diffusivities. In 



Rl/L 


Figure 1. Different implementations of a purely random magnetic field 
(with s = 5/3) lead to consistent values of the isotropic diffusivity for 
Rh/lc <K 1 (with Ic ~ O.IL). Circles and upward triangles are for the 
continuous model, and downward triangles are for the discrete model. For 
Rl > Ic, the diffusivity scales as R^, but the results from the different 
magnetic field models diverge with increasing TJl for the reasons discussed 
in Appendix A. For s = 3/2, similar trends are observed. 


the context of sub-grid modelling of cosmic ray diffusion, the scale 
L is close to the smallest resolved scale in an MHD simulation, 
and the effects of magnetic fields at larger scales, approximated 
by Equation (2), should be faithfully reproduced by the cosmic-ray 
propagation equation with only the resolved magnetic field, using 
the diffusion tensor Kij. 

The continuum and discrete magnetic field models are com¬ 
pared in Appendix A where we present the running diffusion coef¬ 
ficient obtained for a various values of Rh/b and numerical reso¬ 
lutions. The range of variation in the values of the diffusivity, ob¬ 
tained under widely varying numerical resolutions (512^-1536®) 
in model (12), and number of modes (N = 150-1024) in model 
(9), decreases from less than 20 per cent for Ri^/lz = 5 to about 5 
per cent for Rl/Ic < 1. Based on these tests, we initially decided 
to use N = 150 and fcmax/fco = 96 with the continuous model 
(9) and the resolution 512® with fcmax/fco = 128 with the discrete 
model ( 12 ). 


3 DIFFUSION IN A PURELY RANDOM MAGNETIC 
FIELD 

Figure 1 shows the asymptotic (in time) isotropic diffusion coeffi¬ 
cient kq as a function of Rh/L for s = 5/3 and r; = 1 (zero mean 
field) obtained with various magnetic field models and resolutions. 
A linear scaling at 7 ?l <C Ic and a quadratic scaling at i?L Ic are 
found. The latter scaling is consistent with both theory and earlier 
test particle simulations (see Parizot 2004; Plotnikov et al. 2011, 
and references therein). However, the differences between the dif¬ 
fusivities obtained with different magnetic field models increase 
with RhjL for a. fixed resolution and fcmax/fco. As discussed in 
Appendix A, this can be remedied by changing the range, number 
density and/or location of the magnetic modes in fc-space. 

On the other hand, results at Ri^ <C Ic exhibit an encourag¬ 
ingly weak dependence on the details of the magnetic field model. 
Figure 2 illustrates this regime in finer detail where we have added 


MNRAS 000, 000-000 (0000) 




6 A. P. Snodin, A. Shukurov, G. R. Sarson, P. J. Bushby and L. F. S. Rodrigues 




Rh/L R^/L 


Figure 2. The isotropic diffusion coefficient kq at ~ O.IL for 

77 = 1 and s = 5/3 with the best fit to the most highly-resolved data shown 
with (green) solid line (see the text and Figure 3). The linear dependence 
of kq on extends to very small values of _Rl, provided fcmax/fco is 
large enough. The vertical lines of the corresponding colour indicate = 
27r/fcmax for the various parameter values used. The numerical resolutions, 
mode numbers and magnetic field models used are specified in the legend. 
Deviations from the asymptotic linear dependence of kq on Ri^jL at small 
values of this ratio are sensitive to these factors. 


additional results from a discrete magnetic field model at a reso¬ 
lution 1536® with fcmax three times larger than that with the 512® 
resolution. In the 512® case, deviations from the linear trend shown 
by the solid line line start almost from the extreme right of the range 
of i?L shown. The solid blue vertical line (the second vertical line 
from the right) shows where i?L = 27r/fcmax, i-c. at the scale of 
the smallest structures in the magnetic field; significant deviations 
from the linear trend occur below this marker. With the 1536® grid, 
deviations from the straight line are much weaker (the leftmost ver¬ 
tical line marks 7 ?l = STr/fcmax for this case). With the continuum 
magnetic field model, the deviation from the linear trend is weaker 
than that in the case of a discrete model with a similar value of 
fcmax- These results seem to indicate that the continuation of the 
linear trend down to smaller i?L depends on the presence of mag¬ 
netic modes at or near to the resonant wave number fcros = 27r/i?L- 
Suppression of such modes presumably leads to an increased parti¬ 
cle mean free path parallel to the local magnetic field, resulting in 
an enhanced diffusion coefficient, as observed in Figure 2. Given 
that the continuum model has very few modes at large k, it is per¬ 
haps surprising that deviation from the linear trend occurs at ap¬ 
proximately i?L = 27r/fcmax, whereas the discrete magnetic field 
model, containing many more modes near that scale, deviates from 
the linear trend at a somewhat higher value of i?L- One reason for 
this can be that the high-fc modes in the discrete model are not suf¬ 
ficiently well resolved in fc-space. It may also be the case that the 
numerical construction curtails energy at the smallest scales. Fur¬ 
thermore, at these extremes i?L is close to the grid resolution on 
which the magnetic field is defined, and errors from the trilinear 
interpolation between mesh points can have some effect. These re¬ 
sults suggest that caution is required, and that i?L should not be less 
than about three separations of the mesh on which magnetic field is 
defined. 


Figure 3. The isotropic diffusivity kq at low FJl in a continuous and dis¬ 
crete random magnetic fields with the number of modes N, fcmax/fco and 
the resolution specified in the legend. The best fit of Equation (20) is shown 
with dashed line (also shown solid in Figure 2). Open symbols show the 
corresponding relative deviations from the best fit (the right-hand ordinate 
axis). 

Simulations for still smaller values of i?L/T are more compu¬ 
tationally efficient with the continuous magnetic field model. We 
explored the range down to Rl/L = 10“® using up to = 5000 
modes; the diffusion coefficient remains very close to the solid line 
shown in Figure 2. Casse et al. (2001) found a Bohm scaling of the 
diffusion coefficient, k oc Ri^vq, similar to our results in a limited 
range of Ri^. However, these authors find a different scaling below 
a certain value of i?L (see also Fatuzzo et al. 2010). It is not clear 
why the result found here is different from those; the most obvious 
difference is that those simulations used fcmax/fco = 10^ with only 
N = 200 modes in their continuous magnetic field model. Here 
we confirm the linear scaling of kq with Rh/b down to very small 
values of 7 ?l /Ic with both magnetic field implementations. 

In Figure 3 we show a fit of the form 

« 0.0031 + 0.74 ^ (20) 

vqL L 

to the isotropic diffusivity obtained with fcmax/fco = 256 over 
the range where the variation with Ri^jL is linear, 4 x 10~® < 
Ri^/L < 5 X 10~^ with (c ~ 0.1. The relative deviations from the 
fit are shown in Figure 1; they do not exceed about 15 per cent at 
Rh = Ic (i.e., Rh/L « 0.1), where we might have expected more 
deviation, given the tendency towards quadratic scaling with Rh 
at around Rh ~ L- The fit approximates the numerical results to 
within 3 per cent for all values of Rh, between Rh = (c/2 and the 
smallest values of 7 ?l considered (where already Rh < 27r/fcniax, 
resulting in reduced resonant scattering). 

The fit has a non-zero intercept, 

~ 3 X 10"®t;oA, (21) 

corresponding to a part of the diffusivity independent of the Larmor 
radius (i.e., the particle energy). For vo = c, and L ~ 100 pc for 
the outer turbulent scale, we obtain kq — 3 x 10^® cm^ We 
can compare this result to an estimate based on the particle diffu¬ 
sion along random magnetic lines (lokipii & Parker 1968), known 
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as the field line random walk (FLRW) model. If particles are tied to 
the magnetic lines and maintain a constant pitch angle, the diffusion 
coefficient perpendicular to the magnetic field can be estimated as 
the magnetic field line diffusion coefficient (defined as the spatial 
divergence rate of the magnetic lines - see below) multiplied by the 
particle velocity parallel to the field lines, i.e., k = fiDb/S, where 
Df, is the magnetic field line diffusion coefficient, as defined be¬ 
low, and C/3 is the characteristic particle speed along the field line, 
with the factor 1/3 being appropriate for the isotropic case. Note 
that it is more intuitive to make such a calculation when 77 <C 1 , 
where the diffusion coefficient is essentially measured perpendicu¬ 
lar to the mean magnetic field line; nevertheless, here we will asso¬ 
ciate K with perpendicular diffusion, although there is no preferred 
direction in the isotropic case 77 = 1. Let us now define Df, for 
a general 77 , since it will be useful in subsequent sections where 
77 < 1. Consider random magnetic field line displacements 
as a function of the arc length s, measured along the magnetic line 
(i.e. the change in field line position with s). Then defining 


Db,i = I 


lim 


d((Aa:r)^) 

ds 


( 22 ) 


we take Db = \ Db,i for 77 = 1 and Db ^ \ ^b,i 

when 77 < 1 and Bq || $ 3 . Note that Db has the dimension of 
length. We evaluate Equation (22) numerically by solving field line 
equations for many field lines over several magnetic field reali¬ 
sations, essentially as described by Sonsrettee et al. (2015a). Al¬ 
ternatively, Sonsrettee et al. (2015a) provide an analytical approx¬ 
imation of Db in terms of L or another length scale from the 
magnetic power spectrum. In both cases, for 77 = 1 we find that 
Db « 0.06L. Due to the pitch angle scattering, particles travel 
back and forth along magnetic lines, so that v can be much smaller 
than the instantaneous particle speed vq. Equation (21) is obtained 
for V « 0 . 15770 . If this model is correct, then the parallel scattering 
would have to be non-diffusive, otherwise we would obtain com¬ 
pound subdiffusion (Urch 1977; Kota & lokipii 2000) perpendic¬ 
ular to the magnetic lines. A more likely scenario is that particles 
eventually escape to nearby field lines, so that the FLRW limit is 
invalid. 

To gain further insight into the nature of Equation (20), we 
considered various forms of the magnetic power spectrum M(k). 
When its slope is s = 3/2, rather than s = 5/3, we have 


«:o(7?l) 

vqL 


0.0019 -f 0.76 


Rh 


that is, a similar slope with the intercept different by 30 per cent 
from that in Equation (20). We also considered magnetic spectrum 
of the form ( 8 ) with s = 5/3 and kb = 5ko to obtain 


nojRh) 

vqL 


0.0017-f 0.75 




which, again, is similar to the other two forms but with yet another 
value for the intercept, about a half of that obtained for s = 5/3. 
This suggests that depends on the form of the magnetic 

power spectrum, but that this dependence is not very strong, varia¬ 
tions remaining within a factor of two for reasonable forms of the 
spectrum. 


4 DIFFUSION IN A PARTIALLY ORDERED MAGNETIC 
EIELD 



Figure 4. Asymptotic (in time) diffusion coefficients kj_ and Ky, for var¬ 
ious values of 77 specified in the legend, as functions of TJl for the con¬ 
tinuous magnetic field model with s = 5/3 and fcmax/feo = 200. The 
black solid curve shows the diffusivity in the corresponding purely random 
magnetic field as in Section 3. 


Figure 4 shows the local parallel and perpendicular diffusion 
coefficients Ky and k,± for several values of 77 as a function of Ri^, 
for the magnetic spectrum (6) with s = 5/3. A reduction in the 
pitch angle scattering is expected where i?L < STr/fcmax, that is 
Rh/L < fco/fcmax = 5 X 10“® (the latter equality applies for the 
parameter values used in Figure 4). This is evidenced in the de¬ 
viation from the trend of the diffusion coefficients at the smallest 
values of Ri^/L plotted in the figure. Although this is not easy to 
see in the figure itself, the behaviour is similar to that seen in Fig¬ 
ure 2 when plotted on a linear scale. For a given i?L, the diffusivity 
in a purely random magnetic field, kq obtained at 77 = 1 , separates 
regions in Figure 4 filled with the curves representing Ky in the 
upper part and acx in the lower: k,± < kq < Ky. 

Before considering the range Rl/L 1, our main concern 
here, a few comments on the region of larger values of i?L, where 
Rl/L > Ic/L — 0.1 might be appropriate. Here Ky oc while 
Kx tends to become independent of Ri^/L. Plotnikov et al. (2011) 
explored the range 77 > 0.5 to find that Kx/Ky oc R^^ at large 
Rl, even when 77 is very close to unity (which implies a significant 
reduction of k± from kq for any Bo > 0). Our results reproduce 
this behaviour for 77 = 0.5, but we have not explored large enough 
to observe it at larger values of 77 . However, k± decreases with 
increasing i?L for 77 = 0.1 and 0.5, rather than becoming indepen¬ 
dent of Rl/L. In the rest of this section we consider small Rl- 


4.1 Parallel diffusion 

For a given Rl, we find that Ky scales quadratically with Bo/bo — 
or linearly with (1 — 77)/77 — as shown in Figure 5. The numerical 
results can be represented by kq of Equation (20) with the addition 
of a slowly varying power law function of Rl , of the form uRl ■ We 
fitted the two parameters a and b at small values of Rl to obtain 
the following approximation shown in Figure 5 with dashed lines: 


In this section we consider cosmic ray diffusion in a magnetic field 
with 77 < 1 , where Bo is aligned with the z-axis (z = xz). 


«:||(-Rl, 77 ) ^ ko(J7l) ^ 1 / ^ 
voL voL ^ \ L J 77 


(23) 
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Figure 5. The local parallel dilfusivity Acy plotted versus Bo/i>o = 
\/{l — rj)/r} for two values of Ri^, together with their approximations (see 
the text). Filled (blue) circles are for i?L = /c/4, where the relative eiTor 
in the approximation is less than 10 per cent for the range of Bq /5o shown 
(0.1 ^ rj < 1) (as shown by open circles, with the axis on the right of the 
frame). For the smaller value Rj^jL = 0.005 (filled triangles), the error 
(open triangles) is still smaller at a few per cent. 
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Figure 6 . The relative error in the approximation for K|| given by Equa¬ 
tion (23) to the test particle simulations for various values of rj given in the 
legend. The largest values of TJl shown coiTespond to Rl ~ lc/2. 


We note that the final term has the same scaling as is obtained from 
quasilinear theory (Jokipii 1966). In particular, for the magnetic 
spectrum ( 6 ), and assuming magnetohydrodynamic waves propa¬ 
gating parallel to the mean field, the appropriate quasilinear result 
is (e.g., Berezinskii et al. 1990; Schlickeiser 2002, and references 
therein) 


QL _ 


AvoL (27r) 


II (s-l)(2-s)(4-s) 
which for s = 5/3 gives 


BoV ( Ri^'' 


^QL 

^11 

voL 




(24) 


(25) 
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0.05 
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Figure 7. Perpendicular magnetic field line diffusion coefficients from test 
field line simulations as a function of Bo/bo and a theoretically motivated 
approximation (see text) for s = 5/3. Here Dq = Note that 

the factor ADq/Ic ~ 2.4, obtained here analytically, is not expected to 
depend strongly on s. 


which is in very good agreement with the fitted term above. For 
s = 3/2 and rj = 0.5 a fit to the data yields a similar scaling. 
An important point to note is that the usual definition of i?L in 
Equation (24) (i.e., with respect to the mean field) deviates from 
our definition in the limit 77 —>■ 1. Figure 6 demonstrates the fair 
quality of approximation of the test particle simulations for various 
77 and i?L. For i?L < lc/4 « 0.025Z/, the accuracy is better than 
five per cent, and within 10 per cent at Rl = lc/ 2 . 


4.2 Perpendicular diffusion for i?L -C 


We find that for the local perpendicular diffusivity obtained in our 
simulations, the expression 


K±{RL,ri) 


ko{R^) 

l + x(Bo/6o)= 


(26) 


provides a good approximation when i?L -C Ic- Here x is ^ con¬ 
stant around 2 or 3 that is not very sensitive to the value of i?L and 
KoiRh) is given by Equation (20). We have arrived at this form by 
first comparing the magnetic field line diffusion coefficient Di,, as 
discussed in Section 3, with k±, motivated by previous test particle 
simulations (e.g., Michalek & Ostrowski 1997; Ruffolo et al. 2008) 
that have shown that these two quantities are related. We calculated 
D(, for various values of Bo/ho using test magnetic line simula¬ 
tions and compared this with k±_ from our test particle simulations. 
For a given i?L, our results show approximately oc Db for a 
large range of Bo/bo. Then using Equation (20) we can write 


n±{RL,'n) = Ko{RL)Db{ri)/Do, 


(27) 


where Do = Db\Bo=o- As mentioned in Section 3, Do can be 
approximated analytically in units of Ic (Sonsrettee et al. 2015a). 
We have no general expression for Db, but Sonsrettee et al. (2015b) 
have derived analytical expressions in the limit Bo/bo —>■ 00 and 
shown that a simple rational function which interpolates between 
Do and this limit describes test field line simulations very well. 
Following these authors, for large Bo/bo we have Db —>■ D^o = 
(lc/4)(6o/Ho)^, then we obtain the expression 


Dbiv) 


Do 

1 + Do/ Doo 


Do 

1 + (4L»o/(c)(Bo/6o)2’ 


(28) 


n 
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Figure 9. As Figure 8, but with the approximate form of k± given by Equa¬ 
tion (29). 


where we have introduced a rational function that interpolates be¬ 
tween Do and Doo- As shown in Figure 7, this is in very good 
agreement with our test field line simulations. Finally, substitut¬ 
ing Equation (28) into Equation (27), we obtain Equation (26) with 
X = 4Do/(c. If we take Dq obtained directly from field line sim¬ 
ulations we get X ~ 2.35, which is very close to the theoretical 
value X = 2.4 used in Figure 7. 

We also obtained Do from simulation using s = 3/2; for 
these, X ~ 2.44, which only differs from the s = 5/3 value by 
a few per cent, which is about the same size as the error. Moreover, 
since there are reasonable arguments that Do should scale linearly 
with Zc, we expect that x should be essentially a constant, indepen¬ 
dent of the detailed form of the magnetic spectrum. 

In the following subsection we compare Equation (26) and 
some improvements directly with our test particle simulation data. 


4.3 Perpendicular diffusion for i?L < Zc/2 

As shown in Figure 8, x = 4Do/Zc in Equation (26) provides a 
very good approximation to obtained from test particle sim¬ 
ulations. However modest in magnitude, the relative error varies 
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Figure 10. As Figure 8, but with the kj_ given by Equation (30). 


systematically with JZl, inviting a more accurate approximation. 
Figure 8 suggests an additional power-law factor in i?L with the 
exponent of 0.54—0.72 for a given rj. For rj = 0.2, we fitted the 
numerical results in the range 0 < 77 l < b to obtain 

_ 0.19(i7L/L)°» 

voL ~ l + (4I7o/Zc)[(l-77)/r,] ’ ^ ^ 

where the exponent has a standard error of less than one per cent. 
For rj = 0.5, the numerator becomes 0.20(i?Lwith a simi¬ 
larly small error. Furthermore, the numerator exhibits a dependence 
on rj. The relative difference between the fit at 77 = 0.2 and the nu¬ 
merical results at different values of rj is shown in Figure 9. For 
»7 —>■ 1, the difference varies linearly with Ri^/L, similarly to kq 
in Equation (20) or the related approximations of Section 3. How¬ 
ever, variation with a smaller exponent is apparent when 77 —>■ 0. 
The following heuristic form with a linear transition between these 
regimes appears to be a suitable improvement of Equation (29): 


n^{Ri,,rj) 

voL 


rj-^ + 0.19 (1 — rj) 
voR 


Rl 


1 + 


4Do 1 — ?7 


Ic 


n 


(30) 


The deviations of the numerical results from this form are shown in 
Figure 10. 


5 SUMMARY AND DISCUSSION 

The main outcome of our work is a set of approximate expressions 
for the diffusion tensor that apply to cosmic ray propagation in 
a random magnetic field that has sufficient energy at the Larmor 
radius scale to scatter the particles. The latter condition is not a 
generic constraint but rather a result of the limitations of the model 
(for example, our model magnetic field does not include the prop¬ 
agation of magnetohydrodynamic waves, which should relax this 
constraint in a physical manner). Our results can be represented as 
combinations of existing theoretical ideas and analytical estimates, 
in forms which may have been difficult to envisage without detailed 
computations. 

Our main motivation has been to make a first step towards 
a sub-grid model of cosmic ray propagation in the MHD context, 
that captures the dependence of the cosmic ray diffusion on the 
local physical parameters of the plasma. We envisage that the ap¬ 
proximations to the elements of the diffusion tensor derived here 
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may be used with the values of 77 and L obtained from the mag¬ 
netic field parameters near the smallest resolved scales of an MHD 
simulation. It is presumed that the effects of the magnetic field at 
scales larger than the numerical resolution would be accounted for 
by the solution of equations for cosmic-ray dynamics, e.g., in the 
advection-diffusion approximation. 

The main contribution to the large-scale dynamics of interstel¬ 
lar and intergalactic plasmas comes from cosmic ray particles in the 
GeV energy range whose Larmor radius 7 ?l is by far smaller than 
any other scale in the MHD approximation. Therefore, a sub-grid 
model should be based on a description of the cosmic ray diffu¬ 
sion in the limit Rj^/Sx -C 1, where Sx is the numerical resolution 
length. Most studies of cosmic ray diffusion are in the opposite lim¬ 
iting case, where the Larmor radius is comparable to or larger than 
the magnetic turbulent scale. Furthermore, calculations of the dif¬ 
fusion tensor Kij are usually presented in a reference frame con¬ 
nected with the global mean magnetic field. In our case, such a 
mean magnetic field is undefined since there are no reasons to ex¬ 
pect any scale separation in the magnetic fields at length scales of a 
few parsecs. Therefore, we present our results in terms of the local 
diffusivities defined with respect to the local magnetic field. 

We have found the following approximations to the elements 
of the diffusion cosmic ray tensor valid when the particles’ Lar¬ 
mor radius is much smaller than the scale of magnetic turbulence, 
Rl Ic, with Zc the correlation length; we defer comparison of 
the approximations with previous results till later in this section. 

In a purely random magnetic field, fits to the simulation data in 
Section 3 suggest that the isotropic diffusivity can be approximated 
by 


iio{Rh) 


vqL I ai + 02 


Rh 


(31) 


where vo is the particle speed, L is the outer turbulence scale and 
the parameters ai and 02 depend on the magnetic power spectrum. 
The factor 02 changes more weakly with the spectrum than ai, 
at least for the spectra explored. For a power-law spectrum with 
the exponent s = 5/3, we have obtained ai ~ 3.1 x 10“^ and 
02 ~ 0.74, compared to oi ~ 1.9 x 10“® and 02 ~ 0.76 for 
s = 3/2. These approximations remain accurate to within about 15 
per cent up to Rl = Ic- For relativistic particles (uo = c) in a field 
of 5 fj,G with L = 100 pc, these expressions imply a minimum 
parallel diffusion coefficient of (2-3) x 10^® cm^ (in the limit 
of a vanishing mean field, r/ —>■ 1 ). 

In a partially ordered magnetic field, with the mean strength 
Bo and the root-mean-square random part 60 , the parallel diffu¬ 
sivity of the form (23), approximates our numerical results for 
0 < Bo/bo < 1/3 to within about five per cent at Rl = Zc/4 
(and better at smaller values of Rl)- The perpendicular diffusiv¬ 
ity is approximated by Equation (26) with x = 4Do/Zc, where 
Do is the isotropic magnetic field line diffusion coefficient (related 
to the divergence of magnetic lines) of Equation (22) for 77 = 1, 
controlled by the integral length scale of the magnetic spectrum, 
Do/lc « 0.6. 

Although it is not our primary concern in this paper, we have 
obtained the following results for Rl < Zc/2. The perpendicu¬ 
lar diffusion coefficient exhibits an i?L-dependence which is not 
captured by the term containing aiq in Equation (26). For Rl > 
10“®Zc, might be better approximated by Equation (29). There 
is a further, weaker dependence on 77 , which can be allowed for 
using Equation (30). A still better approximation might have the 
constant x in Equation (26) replaced by a function of Rl and/or 


77 . Nevertheless, all the given expressions provide a very good ap¬ 
proximation over the range of Rl considered. 

To illustrate these results, consider an application to MHD 
simulations of the ISM adopting, as typical parameters, Bq = 

2 TtG, 60 = 5 /iG, L = 100 pc, s = 5/3 and Rl = 3.3 x 10^^ cm. 

Recent extensive simulations of the supernova-driven ISM have the 
numerical resolution of order (5x ~ 1 pc. The random magnetic 
field expected at this scale is of the order of 6 | ~ 60 ((5a;/L) ~ 

1 fiG. Our intention is to estimate the effective diffusion coeffi¬ 
cients at this scale. For this purpose we can adopt in our diffusion 
coefficient expressions, Sx and b\sx. as the outer scale and fluctu¬ 
ation strength, respectively. The 77 -dependent term in the expres¬ 
sion for the parallel local diffusivity (23) is then of order 10~^ 
and is independent of Sx. This is about three times larger than 
the value of iio/{voL)\L=Sx. estimated in Section 3, and we find 
K|| « 1.4 X 10^^ cm^ s“^. From Equation (26) we obtain k± « 

3 X 10^® cm^ s“^, and K||/kx ~ 50. The random walk of numeri¬ 
cally resolved magnetic lines should lead to a more isotropic diffu¬ 
sion tensor and when applying the usual lengthscale L and fluctua¬ 
tion strength bo in our expressions we obtain K||/kx « 1.5, which 
is in reasonable agreement with the expectation of isotropization 
at large scales (Berezinskii et al. 1990). Then the cosmic ray dif¬ 
fusivities approximated here can be usefully applied to such MHD 
simulations with bo obtained using the random magnetic field at the 
smallest numerically resolved scale. 

There are other plausible interpretations of the estimated dif¬ 
fusivities. Cosmic ray diffusion is controlled at the resonant scales, 
and a more relevant value of bo to be used could be that at the 
Larmor radius scale. The random magnetic field extrapolated to 
the Larmor-radius scale, — 5 o(/7l/L/)^^^ ~ 4 x 10“® /tG, 
is very weak indeed. Since magnetic field at scales exceeding Sx 
is fully resolved, its effect on cosmic ray diffusion does not need 
to be included into the sub-grid model of diffusion. Then in place 
of Bo, the appropriate effective mean magnetic field could be the 
root-mean-square field at the scale Sx, bsx — 1 uG- Taking these 
values, the T/l- dependent term in the isotropic diffusivity (31) is 
negligible, whereas the 77 -dependent term in the expression for the 
parallel local diffusivity (23) is about 30. Using Equation (30), we 
then obtain Ky/Kx — 10^. In this scenario, the cosmic ray diffu¬ 
sion tensor is highly anisotropic at small scales and cosmic rays 
are partially isotropized when propagating in the random magnetic 
field at scales exceeding 5® ~ 1 pc. Further development and im¬ 
plementation of a sub-grid model for cosmic-ray diffusion will be 
discussed elsewhere. 

The random magnetic field resolved in MHD simulations of 
the ISM would enhance cosmic ray diffusion through the random 
walk of magnetic lines, modelled faithfully by solving simultane¬ 
ously the MHD and cosmic-ray equations with the diffusivity ten¬ 
sor dependent on the magnetic field direction. This process is ex¬ 
pected to be more complicated than usually assumed, and perhaps 
very different from the predictions of existing models of cosmic ray 
propagation in random magnetic fields, because a significant part 
of the random magnetic field in the ISM is produced by compres¬ 
sion in the shocks driven by supernova remnants, and by fluctuation 
dynamo action. Both processes produce random magnetic fields 
represented by widely separated, intense structures against a dis¬ 
tributed background; such a random field is strongly intermittent. 
The structures associated with the magnetic intermittency are mag¬ 
netic sheets and filaments, resulting from shock-wave compres¬ 
sion and the fluctuation dynamo (Zeldovich et al. 1990; Beck et al. 
1996; Wilkin et al. 2007), respectively. Interstellar random mag¬ 
netic fields are expected to have a significant non-Gaussian part. 
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and the random-phase approximation, as employed here and in sev¬ 
eral other works, does not apply to such structures. The magnetic 
structures can affect the cosmic ray transport significantly by in¬ 
troducing Levy flights of the particles, which would be trapped for 
a relatively long time within the structures and move more freely 
over larger distances between them. This picture also applies to 
intergalactic space where the fluctuation dynamo may be active 
(e.g., Akahori & Ryu 2011). The fractional volume of magnetic 
filaments produced in the ISM by the fluctuation dynamo is es¬ 
timated to be of order 0.01-0.1 (Sect. 7.4.2 in Shukurov 2007). 
Fletcher & Shukurov (2007) estimate the mean separation between 
the supernova shocks in the ISM to be of the order of 10 pc and 
argue that polarization observations of the nearby ISM at low ra¬ 
dio frequencies are consistent with this estimate. It is noteworthy 
that the localised magnetic structures have scales by far larger than 
either the relevant Larmor radii or the resolution of MHD simu¬ 
lations of the ISM. Therefore, magnetic field models that are free 
from intermittency, similar to those used here, are quite appropriate 
to model cosmic ray diffusion at those small scales, but not for the 
global diffusion. The exploration of these effects is one of the goals 
of our ongoing simulations of the ISM. 

The expressions for the components of the cosmic-ray diffu¬ 
sion tensor obtained here are broadly consistent with the existing 
theoretical and computational results. The isotropic diffusion co¬ 
efficient given by Equation (31) is the sum of an energy indepen¬ 
dent (for relativistic particles) and a Bohm-like components. While 
previous test particle simulations have shown a Bohm-like scaling 
(e.g., Casse et al. 2001; Fatuzzo et al. 2010; Beresnyak et al. 2011), 
it has been suggested that this scaling only applies down to a certain 
value of 77 l much smaller than L- Here we have found, more pre¬ 
cisely, that the scaling continues down to i?L ~ STr/fcmax, which 
seems reasonable. It is not immediately clear why the previous 
simulations are inconsistent with this result or how to obtain the¬ 
oretically the energy-independent part of the diffusion coefficient. 
Some previous works (e.g.. Globus et al. 2008; Harari et al. 2014) 
have fitted isotropic diffusion results by including the theoretically 
motivated quasilinear E^~‘‘ term, equivalent to the term in 
our parallel diffusion expression. Equation (23). Despite this dif¬ 
ference, Equation (31) varies less than 15 per cent from the expres¬ 
sion of Harari et al. (2014) forL/50 < Rl < L/2, when s = 5/3. 
The parallel diffusion coefficient given by Equation (23) appears to 
be consistent with a quasilinear scaling at low 77 , while the sim¬ 
ple addition of the isotropic diffusion coefficient compensates very 
well for the breakdown of this scaling at —>■ 1. One should 
note that the quasilinear result of Equation (24) is based on an 
MHD wave propagation model in the limit of zero-frequency waves 
(Berezinskii et al. 1990; Aloisio & Berezinsky 2004), whereas a 
calculation more consistent with our underlying magnetic field 
model, i.e., magnetostatic turbulence, usually yields an infinite dif¬ 
fusion coefficient (Fisketal. 1974). [This is a problem with the 
original quasilinear theory, relating to scattering at pitch angles of 
6 = 7 r/ 2 . The problem might be overcome in several ways; see 
e.g. Cesarsky & Kulsrud (1973), Tautz et al. (2008)]. One explana¬ 
tion for this apparent good agreement with a slightly inconsistent 
theory might be that the effect of propagating MHD waves is not 
very important when considering asymptotic diffusion coefficients 
(Fatuzzo et al. 2010). And finally, the local perpendicular diffusiv- 
ity K,± scales almost perfectly with the magnetic field line diffusion 
coefficient. 

Regarding test particle simulations, we have shown that dis¬ 
crete and continuous magnetic field constructions can yield essen¬ 
tially identical results over a broad range of particle energies. For 


magnetic fields on a discrete mesh, our results suggest that caution 
is required when i?L is less than a few mesh points or is a signifi¬ 
cant fraction of L. 
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APPENDIX A: COMPARISON OF MAGNETIC FIELD 
MODELS AND RESOLUTION TESTS 

First, for reference. Figure Al compai'es the time evolution of acq given by 
Equation (17) at various values of i?L- For ^c/8 ^ -Rl ^ 2/c, the subdiffu- 
sive regime is rather short (about ten time units) and the diffusivity quickly 
reaches an asymptotic value. For Rl ^ both the diffusion coefficients 
and the time required to reach the diffusive propagation can be very large. 
With open circles, for some values of Rl, we have also plotted in Figure Al 
the usual asymptotic diffusion coefficient as given by Equation (18), evalu¬ 
ated after various times. It is clear that it is consistent with the running dif¬ 
fusion coefficient at large times. However, the lainning diffusion coefficient 
converges at much earlier time, so is therefore more appropriate for prac¬ 
tical evaluation. This is because the running diffusion coefficient depends 
on the local (in time) change in the particle mean squared displacement, 
whereas Equation (18) is a global average and must be evaluated at lai'ge 
times in order to forget the initial evolution. 

To understand the advantages and limitations of the two models of 



Figure Al. Time evolution of the (isotropic) running diffusion coefficient 
for vanishing mean field, Bq = 0, for various values of Rl as specified 
in the legend. The continuous magnetic field model has been used with 
N = 512 and femax/^o = 256. In the case of Rl = 4Zc,/c/2 and 
Zc/64, open circles show the evaluation of the usual asymptotic diffusion 
coefficient, i.e. Equation (18), after time t/to. We see that this is consistent 
with the nanning diffusion coefficient if evaluated at sufficiently lai'ge time. 
However, for practical evaluation of the diffusion coefficients, clearly the 
running diffusion coefficient is preferable, since it converges at much earlier 
time. 


a random magnetic field and to estimate the numerical resolution required 
to achieve a desired accuracy, we compare the diffusion coefficients ob¬ 
tained with the two models. One might expect that the main deficiency of 
the magnetic field constructions will be an insufficient number of modes 
near the scale kres = 27r/RL, which would likely result in lower than ex¬ 
pected levels of scattering. Thus, we consider a few values of Rl and a 
purely random magnetic field, 77 = 1. Both theory and eaiJier numerical 
results suggest different regimes of diffusion for the cases Rl Ic and 
Rl ^ Zc, so we consider the cases Rl = 5Zc, Ic and Zc/5, where Zc is ob¬ 
tained from Equation (7). Figure A2 shows the running isotropic diffusion 
coefficient kq = (1/3) for 5 = 5/3. We focus on the time inter¬ 

vals where differences between the two models are the strongest, so that the 
relative differences between the various lines are far larger in Panel (a) than 
in (b). Note that here it is necessary to have very good particle statistics so 
that we can attribute differences in the diffusion coefficients to changes in 
the magnetic field resolutions, rather than the statistical error in the evalua¬ 
tion of a diffusion coefficient; we estimate the maximum relative standard 
error in the diffusion coefficients to be less than 1 per cent in all cases. 

For Rl = 5 Zc (Figure A2a), the cosmic rays become diffusive at 
t ^ 30Zo. In this regime, convergence of the diffusion coefficients in a 
continuous magnetic field seems to require a large number of modes N 
per unit of kmax/ko when compared with the other panels at lower Rl- 
On the other hand, large values of kmax/ko may not be required because 
the particles are not dominated by magnetic fields at small scales for a 
sufficiently steep magnetic spectrum. Plotnikov et al. (2011) conclude that 
kmax/ko = 10 is adequate for large Rl- For fcmax/fco = 96, the choices 
N = 150 and 600 lead to a difference of more than 10 per cent in the diffu¬ 
sivity. For N = 512 and kmax/ko = 256, we might expect the diffusion 
coefficient to be smaller than with N = 600 since the correlation length 
is smaller (according to the restricted k evaluation of Equation 7). How¬ 
ever, the mode density seems to be too low to realise this with N = 512, 
but taking N = 1024 is sufficient to give the expected result. For a given 
kmax/ko, we find that the diffusion coefficients always decrease with in¬ 
creasing N. Whereas the continuum model concentrates modes near ko. 
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the discrete model has very few modes at low k, which is perhaps why we 
see a large discrepancy between the 512^ resolution discrete model and the 
continuous ones in Panel (a). By increasing both the mesh separation corre¬ 
sponding to wave number ko and the overall mesh size (denoted by a star 
in the figure legend), it appeal's that this discrepancy can be eliminated. We 
tried increasing the resolution alone, up to 1600^, but this had little effect. 

In Panel (b), for i?L = ^c, too few modes in the continuous model 
of Equation (9) appear to cause a problem. Overall though, there is less 
vai'iation in the asymptotic n between various resolutions than in Panel (a). 
With the resolution 1024^ in the discrete magnetic field model and N > 
150 in the continuous model, the asymptotic diffusivities agree within about 
2 percent over the whole time evolution (although the N = 600 case turns 
out to be somewhat special, leading to a lower value of k than for the other 
values of N). It is also quite satisfactory that significant variation in the 
resolution of the discrete model does not affect the result much. 

The difference between results obtained under the whole range of 
magnetic field implementations is even smaller in Panel (c), for = ^c/5. 
The diffusivities in the three continuous models coincide almost perfectly, 
which suggests we have sufficient resonant scattering. It is perhaps surpris¬ 
ing that although the discrete model contains many more modes near feres 
than the continuous one, it seems to require that femax be several times 
larger than feres for convergence, whereas the continuous one is well con¬ 
verged at femax ~ feres - This may be due to the discrete nature of the model 
and is discussed further in Section 3. 

The results here suggest that the two models will be in almost perfect 
agreement in the limit of high resolution and when feo is sufficiently well 
resolved in the discrete model. Clearly the discrete model is more limited 
in the range of i?L that can be considered and presumably the limitations 
mentioned above will apply to test particle simulations in more realistic 
magnetic fields. For the range of i?L studied, we see no obvious effects 
due to the periodic magnetic field. For continuous magnetic field models, 
some previous works have determined reasonable values of N for paitic- 
ulai' cases, which has typically been quantified in terms of the number of 
modes A required per decade of femax/feO’ he. N = Alog 2 Q[femax/feol- 
Examples of reported reasonable values of A are 25 (Fatuzzo et al. 2010), 
100 (Parizot 2004), and 200-300 (Plotnikov et al. 2011). Considering the 
results of figure A2, A = 25 might be reasonable for i?L = /c/5. However, 
for i?L = 5/c this would give W 50 for femax/feo = 95, which would 
appear to be unreasonably small. Our results for N = 512, femax/feo = 
256 seem reasonable over all panels, and this gives A ^ 213, which is of 
similar size to those previously reported. The point to note here is that a 
given reasonable value of A is not universally applicable and will depend 
on both i?L and the desired range femax/feo- 

Admittedly we have not performed such a detailed study of the con¬ 
vergence for 77 < 1. However, for some values of 77 , and at the extremes of 
i?L, we have observed similar asymptotic trends to those discussed above 
and in Section 3, so we suspect that variation of 77 is not a very important 
consideration, except perhaps in the limit 77 —)• 0 . 
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Figure A2. A comparison of the running diffusion coefficient obtained for 
the continuous and discrete magnetic field models under vaiious numerical 
resolutions in fe-space, specified in the legends, for s = 5/3 and a purely 
random magnetic field, Bq = 0, at various values of Ri^\ (a) i?L = 5/c; 
(b) i?L = /c and (c) i?L = /c/5. For the discrete model (12), feo is taken to 
be equal to two mesh separations in wave number space (see page 4), apart 
from the two cases denoted with a stai' in Panel (a), for which feo is equal to 
four mesh separations. 
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